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The  general  single -period  optimal  portfolio  selection  problem  Is 
the  following:  An.  Investor  wishes  to  Invest  his  wealth  In  certain  risky 
assets,  each  of  which  has  a constant  scale  of  return  that  Is  a random 
variable.  He  could  also  borrow  or  lend  and  the  Interest  rates  for  borrowing 
and  lending  are  assumed  to  be  the  same  and  fixed.  The  latter  Is  referred 
to  as  the  risk-free  (or  safe)  asset.  The  objective  of  the  Investor  is  to 
maximize  his  expected  utility  of  wealth  subject  to  his  budget  constraint 
and  certain  federal  or  other  personal  restrictions. 

In  Lintner  [ 13  ] and  Ziemba  et  al.  [ 18  1,  It  was  shown  that  the 
problem  can  be  solved  by  a two-stage  procedure,  provided  that  the  investor 
has  a concave  utility  function  and  that  the  asset  returns  have  a multinomial 
distribution.  In  stage  1,  one  solves  a certain  fractional  program  having  as 
variables  the  proportions  invested  in  the  risky  assets.  In  stage  2,  one 
uses  the  optimal  solution  obtained  in  stage  1 to  solve  a certain  stochastic 
program  having  one  single  variable  which  represents  the  proportion  invested 
in  the  safe  asset. 

In  a recent  paper,  by  making  certain  special  assumptions  about  the 
covariances  of  the  risky  assets,  Elton,  Gruber  and  Padberg  [ 7 ] have 
derived  some  procedures  for  the  solution  of  the  fractional  program. 
Unfortunately,  their  derivation  was  based  on  a wrong  set  of  necessary  and 
sufficient  conditions  for  the  program. 

Our  objective  in  this  paper  is  to  develop  an  efficient  method  for  the 
numerical  solution  of  the  fractional  program  arising  in  the  first  stage 
the  portfolio  problem.  The  second  stage  is  often  very  easy  to  solve. 

See  [ 18  ].  / 

/ 

The  organization  of  the  paper  is  the  following:  We  first  show / Gy 
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Chat  the  fractional  program  can  be  reformulated  as  a certain  equivalent 
linear  complementarity  problem  which  forms  the  Kuhn -Tucker  necessary  and 
sufficient  conditions  of  optimality  for  a certain  (strictly)  convex 

‘ 

quadratic  program.  Then,  by  establishing  a theorem  which  shows  how  the 
complementarity  problem  can  be  effectively  solved  by  a simplified  version 
of  the  parametric  principal  pivoting  algorithm  as  described  in  Pang  [ 14  ], 
we  derive  an  efficient  algorithm  for  solving  the  problem.  We  shall  also 
outline  how  the  proposed  algorithm  can  be  profitably  applied  to  a specific 
model  with  upper  bounds.  Finally  we  shall  report  some  computational 
results  including  a brief  conqparison  of  our  proposed  algorithm  and  Lemke's 
(see  Lemke  [12  ])• 


T.  THE  LINEAR  COMPLEMENTARITY  EQUIVALENCE 

As  the  case  where  short  selling  is  allowed  is  computationally  much 
easier  to  handle,  this  paper  treats  only  the  case  where  short  selling  is 
prohibited. 

The  fractional  program  under  consideration  can  then  be  stated  as: 
maximize  ^Tx'  / (x,T  Vx')^  (1) 

subject  to  eTx'  - 1 , Cx'  < d and  x'  > 0 . 


Here  V is  the  n x n sysnetric  covariance  matrix  of  the  (stochastic)  returns 
of  the  risky  assets  and  is  assumed  to  be  nonsingular  (or  equivalently, 
positive  definite),  p is  the  n-vector  of  expected  asset  returns  in  excess 
of  the  risk-free  return,  x1  is  the  n-vector  of  proportions  of  wealth 
Invested  in  each  asset,  e is  the  vector  of  1's  and  n is  the  total  number 
of  risky  assets  under  consideration.  The  matrix  C is  m by  n. 


-a 
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The  theorem  below  shows  that  the  program  is  equivalent  to  a linear 
complementarity  problem.  It  is  an  extension  of  a result  established  in 
l 18  ] for  the  special  ease  where  the  constraints  Cx'  < d are  not  present. 

Theorem  1.  Suppose  that  problem  (1)  is  feasible  and  that  i*Tx'  >0  for 
some  feasible  x' . Then  problem  (1)  Is  equivalent  to  the  linear 


complementarity  problem 


u ■ - n + Vx  + (CT  - edT)y  >0  x > 0 


v ■ - (C  - de  )x 


>o  y > o 


T T 
u x * v y • 0 


in  the  sense  that  there  is  a one-to-one  correspondence  between  the 
optimal  solutions  to  (1)  and  the  complementary  solutions  to  (2). 

T 

Remark  1 . The  assumption  ()*  x'  > 0)  has  the  interpretation  that  the 
total  expected  return  is  positive  for  some  feasible  portfolios  x' . 

Remark  2.  The  feasibility  of  (1)  is  important  for  the  equivalence  to 
hold.  In  fact,  simple  examples  can  be  constructed  so  that  (2)  has  a 
complementary  solution  but  (1)  is  infeasible. 


Proof  of  the  theorem.  Since  the  objective  function  in  (1)  is 
homogeneous  in  x',  the  problem  is  equivalent  to 


maximize  |ATx"  / (x,,^Vx")^ 


subject  to  Cx"  < dex"  , x"  > 0 and  x"  i*  0 
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under  the  transformation  x'  - x"/  eTx"  . We  show  that  this  latter  program 
is  equivalent  to  the  convex  quadratic  program  below 

minimise  Jx^x  - ^Tx  subject  Cx  < deTx  and  x > 0 . (2') 

Indeed,  let  x"  be  an  optimal  solution  to  (1*).  We  observe  that 
HTx"  > 0 so  that  the  vector  x - ( p,Tx"  / x,,7Vx")x"  is  feasible  for  (2 ' ) . 
Moreover,  ve  have 

ixTVx  - uTx  - - OiTx") 2 / x"^Vx"  < 0 . 

T T 

Let  y be  a feasible  solution  to  (21)  with  Jy  Vy  - p,  y < 0 . Then  y ^ 0 

is  feasible  for  (t')  and  we  have  p,Ty  > Jy^Vy  > 0 . Thus 

JyV  - ^Ty  > - i(nTy)2  / yTvy  > - i(pTx")2  /x"V"  - - pTi  . 

Therefore  x is  optimal  for  (2').  Conversely,  let  x be  an  optimal  solution 
to  (2').  We  claim  ^Tx  > 0.  Indeed,  if  nTx  < 0 and  if  x is  a feasible 
solution  for  (2')  such  that  p.Tx  > 0 , then  since  the  vector  y * (p,Tx  / x\x)x 
is  also  feasible  for  (2'),  we  have 

0 < Jx^Vx  - k*Tx  < Jy^Vy  - (iTy  • - J(^Tx)2  / x^Vx  < 0 

T- 

which  is  impossible.  Therefore  pb  x > 0 . Similarly,  we  may  deduce  that 
^Tx  - x^Vx  . Now  let  x"  be  a feasible  solution  to  (V)  with  tiTx"  > 0. 

Then  the  vector  x ■ 0iTx"  / x,,TVx")x"  is  feasible  to  (2').  Hence,  we  have 

- i(^Tx)2  / x^Vx  - Jx’vx  - pTx  < fcx’Vx  - |iTx  - - £(p,Tx")2  / x,,TVx" 

which  implies 

l*Tx  / (x'^Vx)*  > nTx"  / (x',TVx")^  . 
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Therefore  x is  optimal  for  (1').  The  proof  can  now  be  completed  by  noting 
that  the  problem  (2)  is  precisely  the  set  of  Kuhn-Tucker  necessary  and 
sufficient  conditions  of  optimality  for  the  quadratic  program  (2'). 


2.  THE  PARAMETRIC  APPROACH 

I 

Among  the  various  algorithms  which  can  be  used  to  solve  the  linear 
complementarity  problem  (2)  is  the  parametric  version  of  Graves' 
principal  pivoting  algorithm  (see  Graves [ 9 ] for  the  original  non-parametrlc 

version  and  Cottle  [ 2 ] for  the  parametric  version).  A typical  pivot 

step  of  the  algorithm  is  outlined  as  follows:  Given  a parametric  linear 
complementarity  problem  in  the  canonical  form 

w - r + X s + Mz 

where  X is  currently  of  some  positive  value  £ such  that  r + Xs  > 0 and 
rk  + ^Sk  * ® f°r  8ome  we  f*r8t  Pivot  \ into  the  basis  and  wfc  out  of 
the  basis  if  m^  9*  0 (known  as  a 1 x 1 diagonal  pivot) . Else  we  increase  z^ 
to  a value  until  some  r^  + Xs^  becomes  zero,  in  which  case  we  pivot 
on  m^  and  m^  (a  2 x 2 block  pivot) . 

It  is  a rather  well-known  fact  that  if  the  diagonal  pivot  entry  is 
positive  in  each  step  and  if  each  pivot  is  nondegenerate  (see  [9  ] for 

the  handling  of  degeneracy) , the  algorithm  always  terminates  with  a solution 
to  the  linear  complementarity  problem 

1 

T 

w ■ r + Mz  > 0 , z > 0 and  z w ■ 0 

in  a finite  number  of  steps  (when  X reaches  zero) . A sufficient  condition 


r 1 
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for  the  diagonal  pivot  entries  to  be  positive  is  that  the  original  matrix  M 
should  have  all  principal  minors  positive. 

The  matrix  M in  problem  (2)  which  is  given  by 


T 

where  A - C - de  , certainly  does  not  satisfy  this  sufficient  condition. 
Nevertheless,  as  the  next  theorem  shows,  the  same  assertion  about  the 
diagonal  pivot  entries  remains  valid  for  the  problem.  The  proof  of  the 
theorem  is  given  in  the  Appendix. 

Theorem  2.  Consider  the  solution  of  the  linear  complementarity  problem 
(2)  by  the  parametric  algorithm  described  above,  where  the  parametric  vector  s 
is  chosen  as  s * . Then  the  diagonal  pivot  entry  in  each  pivot  step  is 

positive.  In  other  words,  the  problem  can  be  solved  by  performing  the  1 xl 
diagonal  pivots  exclusively. 

Remark.  The  positive  definiteness  of  V is  important  for  the  theorem  to 
hold.  In  fact  Graves  [ 9 ] showed  that  all  the  pivots  are  2x2  if  V is 

the  zero  matrix. 

As  explained  in  [ 14  ] , the  only  information  one  needs  to  have  in  order 
to  execute  the  1x1  diagonal  pivots  consists  of  (i)  the  index  set  of  the 
currently  basic  z -variables,  (ii)  the  current  constant  (r  - ) column  and 
(iii)  the  current  parametric  (s  - ) column.  The  update  of  the  matrix  M is 
entirely  unnecessary. 


1 


Applying  this  idea  to  the  problem  (2) , we  may  formulate  the  algorithm 
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below.  (See  the  Appendix  for  the  explanation  of  the  notations  and 
recall  that  A ■ C - deT  with  C being  in  by  n.) 

Algorithm  for  solving  problem  (2) . 

Step  0 (Initialization)  Set  0 - 6 - 0,  or  - {l n}  and  y » {l m}. 

I M ** 

Step  1 (Main  computation)  Solve  the  system  of  linear  equations  for  / Hg 


Step  2 (Ratio  test)  Determine 

\ » max  { max  :ei>0  i €«}  , max  : <0,  j €p  } , 

max  { b|  / b^  : b^  > i 6 y } » { b j ^ bj  : b j < j € ft  } } . 

If  \ < 0,  terminate  with  the  solution 


Otherwise,  let  k be  a maximizing  index  and  continue. 


8 

Step  3 (Updating  the  index  sets)  Let 

'8  U (k)  if  ratio  occurs  at  first  maximum 

8 ■ < P \{k  } if  ratio  occurs  at  second  maximum 

otherwise 

and  a - {t ,. . . ,n  ysjj  . Update  y and  & in  a similar  fashion.  Go  to  Step  1„ 

We  point  out  that  the  major  computational  effort  required  by  the 
algorithm  la  contained  in  Step  1.  In  general,  this  step  should  best  be 
implemented  by  using  an  adaptive  procedure  (such  as  those  described  in 
Gill  et  al.  I 8 ])  to  take  advantage  of  the  change  of  the  index  set  $U6. 

3.  A MODEL  WITH  UPPER  BOUNDS 

In  I 7 ] , by  making  certain  assumptions  about  the  covariance  matrix  V, 
Elton  et  al.  studied  several  special  cases  of  the  following  model 

maximize  p,Tx*  / (x'^Tx')^  subject  to  eTx'  ■ 1 and  0 < x’  < d (5) 
which  of  course  is  equivalent  to  the  complementarity  problem  (2)  with  C being 

the  Identify  matrix  of  order  n.  We  remark  here  that  the  linear  complementarity 
formulation  used  by  the  authors  of  the  aforementioned  reference  is  not  the 
correct  one.  Our  purpose  in  this  section  is  to  outline  how  Step  1 of  the 
proposed  algorithm  can  be  greatly  simplified  by  taking  advantage  of  the  fact 
that  C is  the  identity  matrix.  In  fact,  the  entire  algorithm  can  be  carried 
out  by  operating  on  the  matrix  V and  the  vectors  p,  and  d only. 

Referring  to  the  notations  of  the  algorithm,  we  state  the  following 
result. 

Theorem  3.  Throughout  the  solution  procedure , 5 c p . 


The  idea  of  the  proof  is  in  fact  quite  simple.  By  making  an  inductive 


assumption  that  the  assertion  is  true  before  a certain  pivot  and  then  by 

deriving  an  explicit  expression  £or  the  current  constant  and  parametric 

columns,  it  can  be  observed  easily  that  the  maximum  ratio  (in  Step  2 of  the 

algorithm)  will  not  occur  at  those  rows  corresponding  to  the  basic  variables 

x.  and  v in  the  current  canonical  tableau  of  the  complementarity  system 
8 or 

(cf.  (7)  in  the  Appendix).  Consequently,  the  next  pivot  will  not  occur  in 
these  rows,  thereby  establishing  the  claim  for  this  (and  thus  all  subsequent) 
pivot  step(s). 

Results  similar  to  this  one  have  appeared  in  the  study  of  quadratic  programs 
with  only  upper  and  lower  bounds  on  the  variables  (see  [14])  and  in  that  of  a 
certain  piecewise  linear  complementarity  problem  (see  Kaneko  [11]).  Roughly 
speaking,  8 and  a consist  of  those  indices  whose  corresponding  variables  x£ 

and  x'  are  at  their  upper  and  lower  bounds  respectively,  and  P\8  those  that  are 
a 

between  the  bounds.  The  theorem  then  says  that  a variable  which  is  at  one 
extreme  can  not  immediately  reach  the  other  extreme  without  attaining  some 
intermediate  values.  Intuitively,  this  is  quite  obvious.  The  theorem  also 
permits  one  to  Ignore  certain  rows  in  the  ratio  test.  This  certainly  is  a 
computational  saving.  For  more  discussion,  see  the  above  two  references. 

Consider  now  the  solution  of  the  system  of  linear  equations  (4a) . 

According  to  Theorem  3 , we  may  write  jl  ■ 6 U T|  with  TJ  H 6 - 0 . Defining 


BCM) 


/ v»  <vT\ 

Us  ° ] 


as  in  the  proof  of  Theorem  2 and  recalling  A ■ I - de  , we  may  write 


- e6d 
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where  I ^ denotes  the  identity  matrix  of  order  6 | . By  defining 


v1i« 

°\ 

(° 

van 

v66 

1‘ 

, Bj  (3  ,6)  - 

0 

0 

“X6 

0 1 

u 

N ° 


e6  0 


\° 


it  is  obvious  that 


B(0,«)  -B(9,8)  +£^3,8)  E2(3,8): 


which  shows  that  B(3,6)  differs  from  B(3,8)  by  a rank- two  matrix.  By 
an  easy  calculation,  we  may  deduce 

.-1 


B(0,6) 


-1 


vm 

0 

\-v 


0 

0 


v_1v 


-i. 


-i 

Birin 


V86  ' V611VT1T1VT18 ' 


Hence,  according  to  the  Sherman -Morrison -Woodbury  formula  (see 
Householder  [10,  p.124]  e.g.)  we  have 

B(3 ,6)”^  -B(0,8)‘1  - I1(3,6)G(3,8)‘1E2(0,8)T 


5,  (3,6)  - 5(3,6)  E,  (3,8) 


- d. 


\ 


V-1e 


CV..-V, V.  V~,le  -e. 
66  6T|  TITJ  T|fi  6 BT|  T|T|  T)  6 


where 


B2(3,6)T  - B2(a,fi)TB(S,6)_1  - 


-Ti  in 


\ 


-d^V.  V 


1 


6 8T1  TTH 


G«,8)  - 


- e^V*^e 

71  TTH  11 


4fa«  - VnnV4« 


irwin  » 


,+4?<v™v-» 


J 


Notice  that  this  last  matrix  6(0,8)  is  2x2.  Hence,  we  obtain 


^71  ’ V " % (471  ’ V 
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°vv  - <“»  • v * vvV  • 

Consequently,  the  solution  of  (4a)  can  be  achieved  by  the  following 


steps: 


1 . Solve  the  system  of  linear  equations  for  (d_  , , O 

ti  ~ ’ ir 

WS  * *n  * V “ <viiad5  • *t\  • V <«•) 

and  compute 


<5, , ;s . - w8t<ij . », . - ^<5, . i, , (6b) 


We  point  out  that  the  major  confutations!  effort  required  to  solve 
(4a)  has  now  been  reduced  to  the  solution  of  (6a)  and  the  confutation  of 
(6b) . This  reduction  is  significant  because  (6)  involves  only  the  matrix  V 
whose  size  is  one  half  that  of  M (cf.  (3)). 


Finally,  we  mention  that  similar  savings  can  be  achieved  in  computing  (4b). 
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4.  COMPUTATIONAL  RESULTS 

In  ordtr  to  numerically  test  the  proposed  parametric  approach,  we  have 
implemented  it  for  solving  some  randomly  generated  problems  of  the  type  (5). 

There  are  two  families  of  problems  being  solved.  Each  of  them  is  characterized 
by  a certain  form  of  the  covariance  matrix.  The  reason  that  we  have  chosen  these 
particular  families  in  the  experimentation  is  because  they  are  among  the  most 
comonly  used  models  in  portfolio  analysis.  In  the  first  family,  the 
matrix  has  the  form 

V - £ + L LT 

where  £ is  an  n x n diagonal  matrix  with  positive  diagonal  entries 
and  L is  an  n x m matrix  with  m much  less  than  n.  This  structure  arises 
from  an  m-index  model  ([14]  and  Sharpe  [17]).  The  case  m - 1 corresponds 
to  the  single-index  model  (Sharpe  [16]).  The  constant  correlation 
coefficient  model  studied  in  Elton,  Gruber  and  Padbexg  [5,  7]  also  gives 
rise  to  a covariance  matrix  having  the  above  structure  with  m - 1 . 

In  the  second  family  of  problems  solved,  the  covariance  matrix  is  given 
in  partitioned  form  V « (V^)  where  for  i,  j - 1....N, 


{ 


Ci/(fV 


rl  + c11fl<«1)T 


t * j 

i - j 


where  C - (c^)  is  an  arbitrary  N x N symmetric  matrix  of  scalars,  f1  is 
an  -vector  and  £4  Is  an  ^ x diagonal  matrix  with  positive  diagonal 


entries.  For  N • 1,  this  structure  reduces  to  that  arising  from  the  single- 
index  model.  In  general,  it  assumes  that  the  risky  assets  are  divided  into 
groups  so  that  members  in  each  group  satisfy  the  assumptions  of  a single- 
index  model.  The  multi -group  model  discussed  in  Elton  and  Gruber  [4]  and  Elton, 
Gruber  and  Padberg  [ 6 ] gives  rise  to  a covariance  suitrix  having  this 
structure. 


In  many  practical  applications,  both  m (the  number  of  indices)  and  N 
(the  nunber  of  groups)  are  fairly  small  compared  to  n (the  number  of  risky 
assets).  Advantage  can  be  taken  of  this  fact  to  further  reduce  the 
computational  effort  (and  in  fact  the  computer  storage  as  well)  required  by 
the  proposed  approach  in  solving  problems  with  these  structures.  To  avoid 
complicated  notations,  we  choose  not  to  present  the  technical  details. 

Two  sets  of  experiments  were  performed  on  a DEC-20  computer  at  the 
computation  center  in  Camegle-Mellon  University.  The  computations  were 
done  in  double  precision  to  reduce  round  off -errors.  The  computer  codes 
were  written  in  FORTRAN  and  the  timings  reported  are  exclusive  of  inputs 
and-  outputs . 

The  first  set  of  experiments  was  concerned  with  the  implementation  of 
the  proposed  method  for  treating  an  m-index  model  and  an  N-group  model.  The 
objective  was  to  test  the  capability  and  efficiency  of  the  method  for  solving 
problems  of  considerably  large  size.  The  data  were  generated  as  follows: 

Each  component  of  the  vector  d was  the  same  and  equal  to  1.75/n.  The  number 
1.75  was  used  as  a control  of  the  total  number  of  pivots  and  the  total  number 
of  variables  at  their  upper  bounds  so  that  these  numbers  would  not  become  too 
small.  For  an  m-lndex  model,  die  diagonal  entries  of  Z were  set  equal  to  2.0. 
The  entries  of  the  matrix  L were  generated  randomly  in  (-1.0, 1.0)  and  the 
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components  of  Che  veccor  p.  in  (0,1).  For  an  N-group  model,  the  components  of 

Che  vector  pi  were  generated  randomly  in  (0,10)  and  thoee  of  f*  la  (-1,1). 

The  diagonal  entrlee  of  S1  were  equal  to  2 ran  + 1 where  ran  waa  a random 

T 

amber  between  0 and  1.  Finally,  the  matrix  C was  equal  to  OG  where  G 
waa  an  N by  N random  matrix  whose  entries  were  random  numbers  between  -1 
and  1.  The  results  are  summarized  in  Tables  1-4  below. 


m 

# of  pivots 

# of  variables 
between  bds. 

# of  variables 
at  upper  bds. 

total  CPU 
time  (in  sac.) 

CPU  time/pivot 
(in  sec.) 

5 

410 

120 

49 

21.1867 

0.0516 

10 

402 

80 

72 

30.3223 

0.0753 

15 

414 

42 

92 

41.8303 

0.1010 

20 

432 

34 

99 

56.261 

0.1303 

25 

418 

27 

100 

67.2163 

0.161 

0.1962 

30 

438 

26 

102 

85.852 

Table  1:  Multiple- index  model  n - 200 


m # of  pivots 

# of  variables 
between  bds. 

# of  variables  total  CPU 
at  upper  bds.  time  (in  sec.) 

CPU  time/pivot 
(in  sec.) 

5 1230 

380 

143 

192.4927 

0.1565 

10  1258 

264 

209 

279.7347 

0.2223 

15  1259 

105 

289 

364.0687 

0.2891 

20  1348 

56 

315 

485.9643 

0.3604 

25  1321 

36 

326 

574.019 

0.4344 

30  1290 

29 

330 

673.5277 

0.522 

Table  2:  Multiple- index  model 

n - 600 
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# members 
in  each  grp. 

# of  pivots 

# variables 
between  bds. 

# variables 
at  upper  bds. 

total  CPU 
time  (in  sec.) 

CPU  time/pivot  j 
(in  sec.) 

20 

197 

79 

15 

2.839 

0.0144 

40 

428 

164 

33 

11.685 

0.0273 

60 

607 

243 

47 

24,418 

0.0402 

80 

805 

321 

71 

42.842 

0.0532 

100 

1026 

410 

81 

68.307 

0.0666 

Table  3:  Multiple -group  model  K 

- 5 

# members 

In  each  grp. 

# of  pivots 

# variables 
between  bds. 

# variables 
at  upper  bds 

total  cm 
. time  (in  sec) 

CHJ  time/pivot 
(In  sec) 

20 

811 

317 

73 

51.799 

0.0639 

40 

1656 

628 

141 

188.314 

0.1137 

60 

2439 

955 

216 

399.642 

0.1228 

80* 

3251 

1233 

300 

721.251 

0.2218 

100 

4163 

1579 

353 

1,088.565 

0.2615 

, Table  4:  Multiple-group  model  N » 20 


The  objective  of  the  second  set  of  experiments  was  to  compare  the  proposed 
method  with  Lemke's  algorithm  for  solving  problems  of  the  above  type.  The 
code  that  we  used  for  the  latter  algorithm  was  called  LCPBIG  and  was  written 
at  the  Systems  Optimization  Laboratory  of  the  Department  of  Operations  Research 
at  Stanford  University.  (The  author  is  graceful  to  Professor  R.  W.  Cottle 
for  making  this  code  available.)  The  data  were  generated  exactly  as  above, 
except  that  the  upper  bounds  were  equal  to  1.35/n.  The  comparison  is 


N 

# members  # of 

in  each  grp.  pivots 

t variables 
between  bds. 

# variables 
at  upper  bds. 

total  CPU  time  (in  sec.) 
proposed  method  Lemke 

2 

5 

23 

7 

3 

.055 

.275 

2 

10 

40 

10 

10 

.128 

1.243 

2 

20 

81 

23 

16 

.417 

8.400 

3 

2 

11 

3 

3 

.024 

.111 

3 

10 

61 

15 

14 

.311 

3.678 

4 

10 

87 

23 

16 

.564 

9.080 

5 

8 

89 

23 

17 

.632 

9.612 

17 


summarized  in  Tables  5 and  6 below. 


n 

m 

# of  pivots  # variables 
between  bds. 

# variables  total  CPU  time  (in  sec.) 

at  upper  bds.  proposed  method  Lemke 

20 

2 

39 

11 

8 

.185 

1.107 

20 

14 

40 

10 

9 

.777 

1.110 

30 

2 

63 

17 

13 

.384 

3.930 

30 

18 

56 

6 

19 

2.025 

3.073 

30 

20 

65 

7 

18 

2.352 

4.027 

40 

2 

89 

27 

13 

.766 

9.877 

40 

10 

88 

14 

21 

1.763 

9.850 

40 

20 

98 

12 

21 

4.067 

11.364 

Table  5:  Comparison  for  an  m- index  model 


Table  6:  Comparison  for  an  N-group  model 
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We  point  out  four  remarks. 

1.  The  numbers  of  pivots  shown  In  the  last  two  tables  are  the  same  for 
both  algorithms. 

2.  In  all  cases,  the  numbers  m and  N are  kept  fairly  small  in  order  to 

be  consistent  with  the  smallest  of  the  ratios  m/n  and  N/n  in  many  practical 
applications  of  these  models. 

3.  The  data  in  all  the  problems  solved  are  extremely  dense.  In  fact,  this 
is  an  essential  reason  why  we  have  compared  the  two  algorithms  only  on 

small  problems.  (The  code  LCPBI6  was  written  for  solving  linear  complementarity 
problems  with  matrix  having  six  thousand  or  less  nonzero  entries.) 

4.  The  proposed  algorithm  is  rather  sensitive  on  the  size  of  m and  N 
(whereas  Lemke's  algorithm  is  not). 

From  these  experiments,  we  may  draw  the  following  two  conclusions: 

1.  In  terms  of  computation  times,  the  proposed  algorithm  is  consistently 
superior  to  Lemke's.  In  most  cases  (when  m/n  or  N/n  is  small)  the  former 

is  several  times  faster  than  the  latter.  The  reader  is  referred  to  Pang  et  al. 
[15]  for  a brief  explanation  based  on  operation  counts  of  the  algorithms. 

2.  The  proposed  algorithm  is  capable  of  solving  large  problems  in  a fairly 
efficient  manner. 

Finally,  we  point  out  that  in  addition  to  the  superiority  in  computation 
times,  the  proposed  algorithm  also  requires  a substantially  less  amount  of 
computer  storage  than  Lemke's  algorithm. 


APPENDIX 


Here  we  establish  the  theorem  below.  Before  doing  so,  we  explain  the 
notations  to  be  used.  Let  A be  an  m by  n matrix.  If  or  and  p are  subsets 
of  {l,...,m]  and  {l,...,n}  respectively,  by  A.  we  mean  the  submatrix  of  A 

«p 

whose  rows  and  columns  are  Indexed  by  a and  0 respectively.  If  j €{l,...,m}, 
we  denote  the  j-th  row  of  A by  Aj.  Similar  notations  are  used  for  vectors. 

Theorem.  Consider  the  solution  of  the  parametric  linear  complementarity 
problem 

0 '0 + xo  +c 

uTx  - vTy  - 0 

where  F€Rnxn  is  symmetric  and  positive  definite,  A.£tH&xn  is  arbitrary, 
e is  the  vector  of  ones,  q is  arbitrary  and  b > 0,  by  the  parametric  version 
of  Graves'  principal  pivoting  algorithm.  Then  all  the  diagonal  pivot  entries 
are  positive.  In  particular,  in  a finite  number  of  steps,  the  algorithm 
will  terminate  with  X * 0. 

We  need  the  following  lemma  whose  proof  is  not  difficult  and  thus  omitted. 

Lemma.  Let  F and  A be  as  in  the  theorem.  Let  P and  6 be  nonempty  subsets 
of  {l,...,n}  and  {l,...,m}  respectively.  Then  the  matrix 


<vT 

0 


) 


is  nonsingular  if  and  only  if  the  matrix  A,  has  full  row  rank.  In  this 


the  matrix  is  Pos^-^vc  definite  and 


where 


is  synmetric  and  positive  semi -definite.  Moreover,  det  B(g,6)  is  positive 


Consider  any  pivot  step  of  the  algorithm.  Let  @ 


and  6 denote,  respectively,  the  sets  of  indices  of  the  currently  basic  x 
and  y variables  before  the  pivot.  With  no  loss  of  generality,  we  may 


assume  that  both  0 and  6 are  nonempty.  Let  B(0,6)  be  the  matrix  defined 
in  the  lemma.  It  is  a rather  well-known  fact  from  the  theory  of  pivotal 


algebra  that  B(0,6)  is  nonsingular.  Then  the  current  canonical  form  of  the 


complementarity  system  can  be  written  in  partitioned  form  as 


‘ r 


21 


where  a and  y are  respectively  the  complements  of  g and  6 in  {l , . . . ,n]  and 


We  divide  the  proof  into  four  cases. 

1.  The  diagonal  pivot  entry  occurs  at  a xQ-row.  In  this  case,  the  entry 

P 

is  HO, 6),,  where  j € g is  such  that  H(g,6).e.  ■ e < 0.  Since  HO, 6)  is 
jj  j P J 

symmetric  and  positive  semi -definite,  the  fact  that  H(g,6).eQ  is  nonzero 

J P 

implies  that  HO»&)jj  is  positive  (see  Cottle  [1]  e.g.). 

2.  The  diagonal  pivot  entry  occurs  at  a y ^ -row.  In  this  case,  the  entry  is 

^A6gFgg^Aflg^  ^ where  Since  the  matrix  1 is 

positive  definite,  the  desired  pivot  entry  is  therefore  positive. 

3.  The  diagonal  pivot  entry  occurs  at  a u^-row.  In  this  case,  the  entry  is 


Fu  ' <rje  % > 


where  J €<*.  It  is  the  Schur  complement  of 


i'^il 

B(g ,0)  in  BO U £ j } , 6)  (see  Cottle  [3]).  This  latter  matrix  BO U { J }»  ®)  is 

nonsingular  by  the  lemma.  Hence  the  desired  pivot  entry  is  positive. 

4.  The  diagonal  pivot  entry  occurs  at  a v^-row.  In  this  case,  the  entry  is 

A.aH(g,6)  (A.a)T  where  j €v  is  such  that  i.aH(P,i)ea  > 0.  If  the  pivot  entry 
JP  JP  JP  P 


22 


V 


were  zero,  then  we  would  have  AjpH(3,fl)  - 0 by  the  symmetry  and  positive 
semi-definiteness  of  H(P,6)(see  [1]  e.g.).  But  this  is  impossible. 

Consequently,  we  conclude  that  no  matter  where  the  next  diagonal  pivot 
entry  is,  it  must  be  positive,  ffiis  completes  the  proof  of  the  theorem. 
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